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SUMMARY  of  FINAL  REPORT 


Our  work  begins  with  a  brief  review  of  important  deterministic  and  probabilistic  phenomenological 
soil  modeling  [1].  Behavioral  ranges  and  associated  problems  are  defined  with  respect  to  random 
excitation,  material  modeling,  free-held  behavior,  and  structure-media  interaction.  Next,  a  number 
of  ideas  are  explored  as  possible  frameworks  for  the  phenomenological  modeling  of  soil  constitutive 
behavior,  accounting  for  uncertainties  [2],  and  a  significant  source  of  experimental  data  especially 
useful  in  our  model  development  is  identified  [7].  Considering  various  possible  approaches,  it  is 
concluded  that  a  Markov  framework  provides  the  best  approach  since: 

1.  The  fundamental  axiom  upon  which  all  Markov  theory  rests  is  the  following:  the  probability 
that  a  system  occupies  some  state  at  „ime  t\  later  than  time  t  depends  on  its  disposition  at 
time  t  and  not  on  any  time  earlier  than  t.  This  is  analogous  to  the  approach  of  (deterministic) 
classical  mechanics  where  the  evolution  of  a  system  is  established  given  its  present  state. 
Thus,  the  Markov  state  transition  matrix  appears  to  be  an  ideal  probabilistic  counterpart  to 
the  transfer  matrix  concept  used  extensively  in  deterministic  modeling,  where  this  transition 
matrix  and  corresponding  probability  distribution  of  states  couples  the  system  evolutionary 
dynamics  with  system  properties. 

2.  The  need  to  know  only  the  most  recent  state  of  the  system  is  attractive  also  from  the  compu¬ 
tational  point  of  view.  It  would  be  possible  to  update  the  state  of  the  system  while  retaining 
data  from  only  one  additional  state.  The  complete  history  is  not  needed  since  it  is  effectively 
incorporated  in  this  most  recent  information.  This  property  also  turns  out  to  be  useful  in 
linking  the  method  to  experimental  data. 

The  central  thread  through  our  efforts  is  the  study  of  the  Markov  transition  probability  matrix. 
The  main  purpose  is  to  explain  different  classes  of  dynamic  behavior  and  constitutive  properties  in 
terms  of  the  mathematical  properties  of  transition  matrices.  These  transition  matrices  are  called 
stochastic  matrices,  and  they  have  many  properties;  a  key  property  from  which  all  others  follow  is 
that  row  sums  equal  1,  that  is,  the  probability  of  being  in  any  of  all  possible  states  must  equal  1  if 
the  set  of  all  possible  states  is  all  inclusive.  The  bulk  of  this  work  is  reported  in  two  reports  [15,16]. 

As  an  aside  to  the  focus  of  our  efforts,  ideas  in  fractal  geometries  as  applied  to  the  problem 
at  hand  drew  our  attention,  and  we  put  a  few  ideas  of  our  own  down  on  paper  for  future  reference 
[17].  Given  the  opportunity  and  time,  we  may  go  back  to  these. 

A  parallel  component  of  our  efforts  centered  about  linking  the  Markov  transition  probability 
matrix  to  experimental  data.  This  is  a  key  to  using  Markov  models,  and,  ironically,  a  topic  that 
does  not  appear  to  draw  much  attention  in  the  applied  literature.  How  does  the  analyst/modeler 
take  experimental  data,  extract  in  a  rational  way  its  essence,  and  build  a  transition  probability 
matrix  that  reproduces  that  data  in  a  probabilistic  sense?  This  was  viewed  as  extremely  important 
given  the  above  studies  regarding  the  theoretical  properties  of  stochastic  matrices  and  given  the 
type  of  data  one  can  expect  [7]  to  utilize.  Our  ideas  are  embodied  in  two  papers  [3,4]  which 
essentially  tell  the  same  story,  but  the  first  for  those  concerned  with  soils,  and  the  second  for  a 
general  audience  that  may  be  interested  in  the  method  for  other  applications.  A  practical  technique 
is  developed  for  the  derivation  of  transition  probability  matrices,  given  data,  and  the  derived  matrix 
is  shown  to  computationally  reproduce  the  data.  We  are  likely  to  use  these  ideas  to  transform  an 
in-house  deterministic  analysis  code  into  one  that  can  account  for  uncertain  soil  media.  A  further 
refinement  of  this  experimental-based  derivation  of  Markov  transition  matrices  is  suggested  for 
more  complicated  behavior.  In  addition,  use  of  such  matrices  as  predictive  tools  would  add  insight 
regarding  their  robustness. 
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TECHNICAL  ABSTRACT 


The  focus  of  this  research  effort  is  the  construction  of  a  mathematical  framework  to  help  us  under¬ 
stand  the  role  of  uncertainties  in  the  constitutive/dynamic  behavior  of  soils.  The  approach  chosen  is 
phenomenological  and  probabilistic.  The  phenomenological/probabilistic  model  is  a  logical  choice 
given  the  nature  of  the  soil  medium  and  the  respective  necessity  for  significant  data.  Soil  behavior 
falls  into  the  category  of  irreversible  systems  and  is  modeled  as  such.  While  the  modeling  of  soil 
behavior  is  our  goal,  we  present  here  some  ideas  on  the  necessary  mathematical  framework  and 
attendant  notes  on  application  by  which  this  is  to  be  accomplished.  It  is  recognized  that  there  are 
many  issues  to  be  addressed  and  resolved  before  the  goal  is  achieved,  and  that  the  first  of  these,  a 
mathematical  framework,  will  have  only  been  touched  by  this  work. 

Markov  theory  is  used  to  provide  a  probabilistic  framework  for  the  modeling  of  constitutive 
uncertainties  in  dynamic  soil  behavior.  This  theory  is  attractive  because  it  is  the  probabilistic 
analog  to  the  classical  physics  approach  to  dynamics,  it  is  easily  recast  into  a  computational  form 
and  can  be  rationally  linked  to  data,  and  because  it  coincides  with  one  concept  of  a  stochastic 
constitutive  model.  Only  the  homogeneous  or  stationary  Markov  model  is  considered  here.  It  is 
expected  that  this  assumption  will  not  severely  limit  the  applicability  of  the  theory. 

In  addition  to  the  Markuv  probabilistic  framework,  a  statistical  framework  is  needed  to  connect 
the  theoretical  model  and  the  data.  This  connection  falls  within  estimation  theory.  While  possible 
approaches  were  maximum  likelihood  estimates,  Bayesian  inference  techniques,  and  Maximum  En¬ 
tropy  ideas,  a  graphical,  data-based  technique  was  adopted  as  providing  the  most  information  with 
a  minimum  of  assumptions.  This  connection  with  data  is  based  on  the  definition  of  the  state  space 
within  which  the  transition  matrix  of  the  Markov  model  operates.  The  mutually  exclusive  and 
complete  definition  of  a  state  space  is,  in  general,  a  very  difficult  task,  given  the  complexity  of  the 
physical  problem  being  addressed,  but  the  approach  chosen  here  appears  reasonable  for  processes 
where  data  is  obtainable.  Thus,  our  work  may  be  roughly  divided  into  a  theoretical  component 
[2,15,16]  and  an  experiment-based  component  [3,4].  That  is,  (i)  given  a  stochastic  matrix,  draw 
conclusions  regarding  the  system  it  represents,  and  (i»)  given  specific  system  information  in  the 
form  of  experimental  data,  translate  the  data  into  the  computational  transition  probability  matrix. 

Some  preliminary  technical  aspects  are  summarized  next:  Our  work  deals  with  stochastic 
matrices  of  order  n,  where  these  are  nxn  matrices  having  each  element  greater  than  or  equal  to 
zero  and  such  that  row  sums  equal  1;  in  symbols: 


A  —  [fl»j]ij=i*  atj  >  0,  Qij  —  1. 

i= i 

Such  matrices  arise  naturally  in  the  study  of  Markov  chains  [9,10]. 

Our  theoretical  approach  is  mainly  geometric  rather  than  algebraic.  Particularly,  we  conceive 
of  a  stochastic  matrix  of  order  n  in  two  different  ways.  Such  a  matrix  A  can  be  regarded  as 
an  operator  uu  Rn,  i.e.,  if  v  is  a  row  vector  in  R",  then  v  -*  v A  defines  a  linear  operator  map 
Rn  -♦  Rn.  However,  by  virtue  of  the  fact  that  A  is  stochastic,  even  more  can  be  raid.  For,  if 
we  allow  A  be  the  simplex  in  Rn  with  vertices  e;(t  =  l,...,n),  where  e,  denotes  a  vector  whose 
iih  component  is  1  and  all  others  are  zero,  then  A  maps  A  into  A.  A  second  way  of  regarding 
the  stochastic  matrix  A  is  to  suppose  it  to  be  a  point  in  Rn  ;  we  identify  the  n 2  elements  of  A 
with  the  n2  components  of  a  point  in  Rn  .  Then  all  the  stochastic  matrices  form  a  closed  bounded 
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polytope  in  Rn  ,  which  we  denote  here  by  E,  a  polytope  being  in  dimension  n  what  polygons  and 
polyhedrons  are  in  dimensions  2  and  3,  respectively.  The  properties  of  A,  which  lead  naturally  to 
the  classification  of  all  stochastic  matrices,  are  intimately  related  to  the  geometrical  description  of 
E. 

The  foremost  property  of  a  stochastic  matrix  is  its  convergence  characteristics. 
Convergence  can  be  defined  with  respect  to  the  Euclidian  norm  in  R"  or  Rn  .  In  the  first  case, 
let  uo  be  an  arbitrary  vector  in  A  C  Rn.  Then  we  can  study  the  existence  of  limm_oouoAm  . 
Alternatively,  regarding  A  as  a  point  of  Rn  ,  we  can  examine  the  existence  of  Iimm_oo  Am.  It  has 
been  demonstrated  that  both  viewpoints  are  closely  related.  An  important  question  addressed  is 
the  relationship  between  the  convergence  properties  and  the  eigenvalues  of  a  stochastic  matrix  A. 
The  main  instrument  for  this  is  the  Frobenius-Perron  Theorem  [18],  for  which  a  proof  is  given  using 
geometrical  methods.  The  main  result  is  that,  given  an  arbitrary  vector  vo  in  A,  limm_oo  voAm 
exists  if  and  only  if  A  has,  for  eigenvalues,  no  roots  of  unity  other  than  1.  However,  to  determine 
more  precisely  the  convergence  properties  of  a  stochastic  matrix,  we  examine  the  existence  of,  and 
characterize  the  properties  of,  limm_00  Am.  If  A  is  an  eigenvalue  of  A,  then  Am  will  be  likewise 
of  Am;  thus,  we  find  that  limm—oc  Am,  if  it  exists,  will  have  for  eigenvalues  only  0  and  1.  It  will 
be  therefore  possible  to  classify  limit  stochastic  matrices  (matrices  of  the  form  limm_00  Am)  via 
the  multiplicity  of  the  eigenvalue  0  or  1.  It  may  seem  that  this  problem  is  complicated  by  the 
fact  that  a  stochastic  matrix  need  not  be  diagonalizable.  (We  may  recall  that  a  diagonalizable 
matrix  A  is  defined  by  the  existence  of  a  non-singular  matrix  V  and  a  diagonal  matrix  A  such  that 
VAV-1  =  A.)  But,  we  show  that  the  eigenvalue  1  always  occurs  in  a  diagonalizable  manner,  while 
0  occurs  likewise  in  limit  stochastic  matrices. 

Also  of  interest  is  the  class  of  nearly-completely  decomposable  stochastic  matrices.  Briefly,  a 
matrix  is  said  to  be  completely  decomposable  if,  by  renumbering  if  necessary,  it  can  be  written  in 
the  form 


A  = 


A2 


AJ 


where  the  A,  are  matrices  which  run  down  the  diagonal  of  A,  all  other  elements  being  0.  Of  course, 
if  A  is  stochastic,  each  A,  will  be  likewise;  hence,  a  study  of  A  simply  reduces  to  the  study  of 
smaller  matrices,  A;,(t  =  l,...,m).  Courtois  [8]  has  studied  stochastic  matrices  that  have  their 
elements  close  to  completely  decomposable  matrices.  Apparently,  these  “nearly-completely  decom¬ 
posable”  matrices  have  many  modeling  applications  using  Markov  processes.  In  connection  with 
such  matrices,  we  prove  the  following  result:  If  B  is  an  arbitrary  limit  matrix,  B  =  limn,—,*,  Am, 
where  A  is  a  stochastic  matrix,  then  there  exists  a  nearly-completely  decomposable  stochastic  ma¬ 
trix  Bj,  such  that  B  =  limm_oo  By,  and  in  fact  Bi  may  be  chosen  arbitrarily  close  to  a  completely 
decomposable  matrix. 

Complementing  the  above  work  on  the  properties  of  stochastic  matrices,  an  experimental-based 
procedure  is  developed  to  derive  Markov  transition  probability  matrices.  For  a  given  experiment 
data  base,  a  state  space  is  defined  and  used  to  derive  transition  probabilities  by  applying  the 
frequency  interpretation  of  probability.  This  approach  can,  in  principle,  be  used  regardless  of 
the  complexity  of  the  behavior,  to  whatever  accuracy  necessary  using  a  more  refined  state  space 
definition.  The  resulting  transition  probability  matrix  is  representative  of  the  behavior  of  the 
medium,  that  is,  the  data  from  which  it  was  derived.  It  is  expected  that  the  matrix  is  reasonably 
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robust  and  representative  for  “similar”  types  of  soils. 


PROBLEM  STATEMENT 


The  focus  of  this  research  effort  is  the  construction  of  a  mathematical  framework  to  help  us  under¬ 
stand  the  role  of  uncertainties  in  the  constitutive/dynamic  behavior  of  soils.  (For  some  background 
on  the  constitutive  modeling  of  soils,  see  [1].) 

A  constitutive  model  or  equation  describes  the  macroscopic  behavior  resulting  from  the  internal 
constitution  of  a  material;  basically,  the  model  characterizes  the  individual  material  and  its  reaction 
to  applied  loads.  A  stochastic  constitutive  equation  therefore  characterizes  the  material  where 
uncertainties  exist  about  its  material  and  geometric  properties,  and  its  reaction  to  applied  loads. 
As  a  simple  case,  consider  an  ideal  elastic  solid  undergoing  axial  strain  resulting  in  axial  stress 
according  to  the  constitutive  relation  a  =  E(  where  E  is  Young’s  modulus  of  elasticity.  If  E  is  a 
random  variable  (or  a  field,  that  is,  a  function  of  position),  then  one  or  both  a  and  e  can  at  best 
be  described  using  probabilistic  descriptors.  If  cr  =  Ec  is  incorporated  into  an  equation  of  motion, 
then  a  stochastic  constitutive  dynamic  model  is  created  by  way  of  the  random  properties  of  E. 

One  may  view  a  generalized  stochastic  constitutive  equation  as  one  which  probabilistically 
defines  how  a  material  undergoes  dynamic  behavior.  A  constitutive  model  must  accurately  describe 
the  experimental  data  used  to  specify  it,  as  well  as  predict  behavior  under  conditions  not  covered 
by  the  original  data.  In  the  development  of  a  constitutive  model,  it  must  be  recognized  that 
experimental  data  is  always  difficult  and  sometimes  impossible  to  obtain.  Therefore,  the  model 
should  contain  as  few  parameters  requiring  evaluation  as  possible.  Furthermore,  the  model  should 
be  a  function  of  the  major  system  variables  so  that  the  underlying  physics  of  the  media/process  is 
adequately  represented,  and  thus  hopefully  understood. 

Of  course,  no  single  constitutive  model  can  represent  any  material  in  all  situations.  Even 
water,  which  is  probably  the  most  studied  and  best  understood  real  material  known  to  man,  is 
never  described  by  a  single  constitutive  law  to  cover  all  situations.  Thus,  whenever  a  constitutive 
model  is  developed,  only  those  features  of  material  behavior  relevant  to  the  physics  at  hand  should 
be  included.  Any  completely  general  formulation,  while  philosophically  pleasing,  will  only  be  of  a 
formal  nature. 

Our  effort  has  centered  about  the  development  of  phenomenological,  as  opposed  to  microstruc- 
tural,  stochastic  constitutive  models  for  soils.  Phenomenological  models  require  data  not  only  to 
determine  parameter  values,  but  also  to  validate  the  model  itself.  Such  model  building  uses  the 
following  procedure  [5]: 

1.  a  generic  form  of  the  model  is  chosen, 

2.  one  examines  families  that  possess  the  general  features  known  about  the  phenomena,  and  picks 
the  class  of  models  that  encompass  the  phenomena  in  a  reasonably  complete  manner, 

3.  evaluate  the  parameters  of  the  model, 

4.  since  the  phenomenological  model  is  probabilistic,  results  will  also  be  probabilistic:  probabili¬ 
ties  of  specific  events,  averages,  sample  function  behavior, 

5.  validate  the  model  using  a  variety  of  data. 

The  phenomenological  model  provides  an  ideal  framework  where  concern  exists  about  the 
correlation  of  real-life  data  with  the  predictive  accuracy  of  the  theoretical  model.  It  also  permits 
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the  analyst  to  retain  the  theoretical  model  while  adopting  the  new  information  that  becomes 
available  with  more  recent  data. 


AN  APPROACH  to  a  COHERENT  FORMULATION 


Markov  theory  has  been  chosen  to  provide  the  framework  for  modeling  constitutive  uncertainties 
in  soil  behavior.  There  are  several  basic  reasons  for  the  adoption  of  this  powerful  theory  in  these 
studies: 

1.  The  fundamental  axiom  upon  which  all  Markov  theory  rests  is  the  following:  the  probability 
that  a  system  occupies  some  state  at  time  tj  later  than  time  t  depends  on  its  disposition  at 
time  t  and  not  on  any  time  earlier  than  t.  This  is  analogous  to  the  approach  of  classical 
(deterministic)  physics  where  the  evolution  of  a  system  is  established  given  its  present  state. 
Governing  differential  equations  are  such  an  example. 

2.  The  need  to  know  only  the  most  recent  state  of  the  system  is  attractive  also  from  the  compu¬ 
tational  point  of  view.  It  would  be  possible  to  update  the  state  of  the  system  while  retaining 
data  from  only  the  previous  state.  The  complete  history  is  not  needed  since  it  is  effectively 
incorporated  in  this  most  recent  information.  This  property  is  also  valuable  in  translating 
experimental  data  into  model  parameters,  as  is  discussed  below. 

3.  Finally,  the  Markov  state  transition  matrix  appears  to  be  an  ideal  mathematical  counterpart 
to  a  concept  of  a  stochastic  constitutive  model,  where  this  transition  matrix  and  correspond¬ 
ing  probability  states  can  be  interpreted  as  a  stochastic  constitutive  model  coupled  with  the 
evolutionary  system  dynamics. 

The  Markov  framework  is  a  conceptual  structure  that  we  use  to  keep  an,  albeit  sophisticated, 
accounting  of  how  uncertainties  in  the  system  or  inputs  propagate  in  space  and  time.  As  a  the¬ 
oretical  construct  it  is  well  understood.  When  it  is  used  to  gain  physical  understanding  about  a 
specific  system,  such  as  a  soil,  it  is  necessary  that  certain  parameters  particular  to  that  system 
be  incorporated  into  its  mathematical  structure.  These  parameters  will  be  representative  of  the 
stochastic  material  and  geometric  characteristics  of  the  system  being  studied.  In  summary,  the 
proposed  model  is  of  two  parts: 

1.  a  probabilistic  theoretical  framework,  and 

2.  a  data-based  parameter  estimation  technique. 

Item  1  reflects  our  understanding  of  how  a  soil  medium  responds  to  its  environment;  we  adopt 
Markov  theory.  Item  2  reflects  the  need  to  connect  the  theory  to  specific  problems.  This  means 
the  utilization  of  data. 

It  is  noted  that  the  behavior  of  a  soil  medium,  under  realistic  conditions,  is  irreversible.  This 
means  that  it  is  not  possible  to  regain  earlier  states  except  in  very  limited  cases.  This  understanding 
has  an  effect  on  the  mathematical  modeling,  and  the  interpretations  of  predictions. 


THE  GEOMETRIC  THEORY  of  STOCHASTIC  MATRICES 


Stochastic  Matrices  as  Operators  on  a  Simplex 

Here,  the  basic  properties  of  stochastic  matrices  conceived  as  linear  operators  acting  on  Rn  are 
defined  and  proved.  Actually,  it  is  seen  that  the  stochastic  matrices  do  even  better:  they  operate 
on  a  convex  compact  simplex  of  dimension  n  -  1,  which  is  contained  in  Rn.  In  fact,  this  simplex  is 
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the  set  of  stochastic  vectors,  i.e.,  the  set  of  probability  distributions  of  a  system  having  n  possible 
states. 

Two  methods  are  involved  in  this  study.  The  first  is  to  identify  the  simplest  types  of  stochastic 
matrices,  determine  their  properties,  and  then  deduce  the  properties  of  the  general  stochastic  matri¬ 
ces  as  built  up  from  these  simplest  units.  (The  simplest  stochastic  matrices  are  the  indecomposable 
ones.)  The  second  method  is  to  relate  the  properties  of  stochastic  matrices  to  their  eigenvalues. 

The  definition  of  decomposability  can  be  motivated  as  follows.  The  numbering  of  the  com¬ 
ponents  of  a  vector  is  completely  arbitrary;  hence,  the  significant  properties  of  a  matrix  which 
represents  an  operator  on  vectors  should  be  defined  so  as  to  be  independent  of  such  numbering. 
Specifically,  let  v  be  a  row  vector  in  R"  and  suppose  v'  is  obtained  from  v  by  some  permutation  of 
the  components.  We  can  write  this  as  v  =  v'T,  where  T  is  an  nxn  matrix  having  only  0  and  1  as  its 
elements;  in  fact,  T  is  doubly  stochastic  and  has  an  inverse,  which  is  also  a  stochastic  matrix.  If  A 
is  a  general  stochastic  matrix  and  w  =  vA,  we  can  write  this  equation  in  terms  of  permuted  vectors 
w',v'  :  w'  =  t/TAT-1.  Thus,  we  see  that  any  significant  property  of  A  should  also  be  shared  by 
TAT-1,  where  T  runs  through  all  possible  permutation  matrices;  rote  that  these  matrices  TAT-1 
will  also  be  stochastic. 

The  definition  of  indecomposability  follows  from  the  above  criterion:  a  stochastic  matrix  A  is 
decomposable  if  there  exists  some  permutation  matrix  T  such  that  TAT-1  has  the  form 

An  A12 

.  0  A22J  ’ 

where  An,  A22  are  ni  and  n2  order  square  matrices,  such  that  nx  >  0,  n2  >  0.  A22  must  necessarily 
be  stochastic,  but  An  is  so  if  and  only  if  Ai2  =  0.  By  means  of  such  permutation  matrices,  it  is 
possible  to  assume  TAT-1  to  be  of  the  form 

Ai 

*  .  > 

.  0  Ar. 

where  A*+i,..-,Ar  are  stochastic  and  indecomposable.  (We  begin  with  the  subscript  k  +  1  since 
there  may  be  k  decomposable  and  stochastic  matrices.)  Indecomposable,  stochastic  matrices  are 
the  easiest  to  study  since  they  turn  out  to  have  a  simple  eigenvalue  equal  to  1.  If  no  other 
eigenvalue  has  absolute  value  1,  then  the  matrix  has  very  simple  asymptotic  properties  which  are 
investigated.  More  generally,  the  properties  of  an  arbitrary  stochastic  matrix  can  be  inferred  from 
its  indecomposable  components. 

The  above  definition  of  indecomposability  was  an  algebraic  one.  There  is  a  geometric  equiva¬ 
lent.  It  was  noted  before  that  a  stochastic  matrix  mapped  a  simplex  of  Rn  into  itself;  we  can  define 
a  stochastic  matrix  by  this  property.  Further,  we  show  that  the  property  of  indecomposability 
translates  into  the  property  that  no  proper  part  of  the  boundary  of  this  simplex  is  mapped  into 
itself.  This  alternative  definition  will  often  turn  out  to  be  more  useful. 

Four  theorems  on  the  eigenvalue  properties  of  stochastic  matrices  have  been  stated  and  proved 

[15]. 

Representation  of  Stochastic  Matrices  as  Convex  Polytopes 

In  the  last  section,  stochastic  matrices  were  viewed  as  operators  acting  on  vectors,  more  particu¬ 
larly  acting  on  stochastic  vectors.  However,  it  is  desirable,  when  considering  the  classification  of 
stochastic  matrices,  to  view  them  also  as  vectors  with  n2  components.  In  fact,  it  is  necessary  to 
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view  them  in  this  way  to  give  meaning  to  statements  such  as:  two  stochastic  matrices  are  close  to 
one  another.  If  these  matrices  are  viewed  as  points  in  R"  ,  the  Euclidian  distance  or  any  topo¬ 
logical  equivalent  enables  us  to  determine  the  distance  between  stochastic  matrices  and  to  define 
continuous  functions  on  the  set  of  stochastic  matrices.  This  then  is  the  viewpoint  taken  in  the 
following  text. 

It  is  then  seen  that  the  stochastic  matrices  form  a  simple  figure  in  R"  ,  namely,  a  convex 
bounded  polytope.  Also,  those  stochastic  matrices  which  are  limits  of  the  form  ltmm_0OAm  can 
be  determined  insofar  as  their  geometrical  status  within  this  polytope.  As  a  final  step,  the  relation 
of  nearly  decomposable  matrices  to  general  stochastic  matrices  is  examined.  Eight  theorems  on 
stochastic  matrices  are  stated  and  proved. 

Consider,  for  example,  the  propagation  of  a  single  wave  front  traveling  in  one  dimension,  the 
real  line,  which  is  discretized  into  intervals  of  equal  length.  Each  interval  is  labeled  by  an  integer 
such  that  consecutive  intervals  bear  consecutive  labels.  Let  u<  denote  the  probability  that  the  wave 
front  lies  in  the  itk  interval  and  is  traveling  to  the  right;  v'{  the  probability  it  lies  in  the  itk  interval 
and  is  traveling  to  the  left.  Let  ve  be  the  row  vector  (. . . ,  tr» _ i ,  ,  t;<+i, . . .)  and  vr  similarly  is 
given  by  (. . . ,  v'i_l,  v[,  v'i+l , . . .).  Let  t;  =  (v/,vr).  Note  that  EjU;  +  E,r'  =  1,  where  all  V{ ,  v\  are 
non-negative.  Given  v  at  time  to,  we  wish  to  calculate  the  transition  matrix  R  such  that  vR  will 
be  the  probability  distribution  of  the  wave  front  at  the  succeeding  time  interval  tj.  Initially,  we 
assume  the  velocity  of  propagation  to  right  or  left  to  be  such  that  the  wave  front  can  travel  the 
length  of  one  interval  in  time  tj  -  to- 

Further,  we  define  that  at  each  boundary  between  intervals  the  wave  front  has  probability  p  of 
being  reflected  and  r  of  being  transmitted;  so  that  p,r  are  non-negative  and  add  to  1.  Suppose  t, 
denotes  the  row  vector  whose  components  are  all  zero  except  the  ith  which  equals  1.  If  v  =  (e;,0)  at 
t  =  t0,  then  v  =  (rei+1,pei)  at  f  =  tj.  In  a  like  manner,  it  is  easy  to  see  that  (0,e.)  is  transformed 
to  (pej,re,_i).  Thus,  we  see  that  R  is  given  by: 


R  = 


rT 

pi 


Pi  ' 

tT~  1  ’ 


where  I,T  are  matrices  indexed  by  ZxZ  :  I jk  =  Sjk,  the  Kronecker  6,  Tj,*  =  &j+i,k,  =  6j-i,k- 

Note  that  T  •  T-1  =  T-1  •  T  =  I,  so  that  the  notation  follows  standard  convention. 

If  we  now  introduce  the  probabilities  {p*},  Pk  -  1>  where  p*  is  the  probability  that  the 

velocity  of  propagation  is  k  intervals  per  unit  time,  we  obtain  the  transition  matrix 


Q  =  £p*R*. 

k=l 


A  number  of  practical  comments  are  due  here: 

1.  the  probability  distribution  of  propagation  velocity  will  be  a  function  of  the  uncertainties 
associated  with  material  properties,  which,  in  turn,  will  depend  on  the  type  and  magnitude  of 
loading;  as  the  material  undergoes  irreversible  deformations,  the  distribution  generally  changes, 

2.  reflection  and  transmission  coefficients  p  and  r  are  measures  of  material  homogeneity  and 
geologic  structure,  and 

3.  we  expect  that  for  larger  k  (extreme  velocities  with  propagation  over  many  elements),  p*  will 
become  exceedingly  small  and  practically  negligible  when  compared  to  the  early  terms  in  the 
series;  thus,  the  infinite  series  can  be  effectively  truncated,  and  in  some  applications,  after  just 
a  relatively  few  terms. 
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We  resume  the  above  development  by  calculating  the  powers  of  R  as  follows:  Let 


then,  clearly  W2  =  I,  and 


WU  =  U“!W, 


where 

Thus, 


R  =  rU  +  pW, 

R2  =  (rU  +  pW)2 

=  r2U2  +  p2I  +  rp(U  +  U"l)W, 

R3  =  r3U3  +  rp2  (2U  +  IT1)  +  [r2p(U2  +  U“2  +  I)  +  p3]W, 


etc. 

It  is  possible  to  modify  the  above  problem  to  obtain  a  finite  order  transition  matrix.  Consider 
now  that  the  boundaries  between  the  -1  and  0  intervals  and  the  n  and  n  4-1  intervals  are  imper- 
miable;  these  boundaries  permit  reflection  only.  The  transition  matrix  R„  for  the  intervals  0  to  n 
is  now  a  stochastic  matrix  of  order  2n  of  the  form: 


Note  that  Rn,  and  therefore  also  Q„  =  w^ere  ^  Pk  are  non-negative,  and 

Pk  -  1,  are  doubly  stochastic.  Thus,  barring  exceptional  cases, 

lim  voQn  = 
k—oo  in 

where  «o  is  an  arbitrary  2n-stochastic  vector,  i.e.,  an  arbitrary  probability  distribution  of  wave 
fronts.  This  means  that  after  a  sufficient  interval  of  time,  the  wave  fronts  are  evenly  (uniformly) 
distributed  in  the  medium,  regardless  of  the  initial  distribution. 

Two  exceptional  cases  are  as  follows: 

1.  r  =  0,  p  =  1, 

2.  t  =  1,  p  =  0. 
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These  cases  can  result  in  Rn  being  a  permutation  matrix.  In  case  1,  Rn  has  n  linearly 
independent  eigenvectors  for  eigenvalue  1,  namely  (e^e*)  for  t  =  1,2, This  means  that  the 
wave  fronts  in  the  ith  interval  at  t  =  0  is  trapped  there  for  all  future  time.  In  case  2,  R*n  is  the 
identity  matrix;  after  2n  time  steps  the  system  resumes  its  initial  distribution  and  is  recurrent. 

In  the  example,  it  is  noted  that  the  finite  order  transition  matrices  were  doubly  stochastic. 
A  consequence  of  this  was  that  regardless  of  the  initial  probability  distribution,  the  limiting  dis¬ 
tribution  was  generally  homogeneous.  Of  course,  this  result  should  come  as  no  surprise  since  the 
transition  matrix  in  each  case  was  derived  on  the  assumption  of  geometric  homogeneity.  Thus,  we 
should  understand  double  stochasticity  of  the  transition  matrix  to  be  a  manifestation 
of  geometrical  homogeneity  of  the  medium  being  modeled.  If  this  assumption  is  dropped, 
we  would  obtain  a  transition  matrix  which  is  not  doubly  stochastic. 

Also,  in  the  expression  for  Rn,  replace  the  t's  and  p's  by  T{  and  p,,  where  i  varies,  but  such 
that  pi  -f  Tj  =  1;  and  so  obtain  a  matrix  which  is  not  doubly  stochastic.  In  fact,  consider  the 
extreme  case  where  all  motion  to  the  left  is  restrained  by  perfectly  reflecting  barriers;  thus,  Rn  has 
the  form 


The  unique  stochastic  eigenvector  of  Rn  for  eigenvalue  1  is  |( en,en );  from  this  we  conclude 
the  expected  result  that  as  t  -*  oo,  all  the  wave  fronts  congregate  in  the  farthest  right  subinterval. 
Note  further  that  since  Rn  maps  (en,0)  to  (0,en)  and  vice  versa,  Rn  is  decomposable  and  has  -1 
for  an  eigenvalue  so  that  R™  does  not  converge  as  m->  oo. 

STOCHASTIC  MATRICES  and  DYNAMIC  PROBLEMS 

n 

Consider  substochastic  matrices  of  order  n:  P  =  (p,j)"J=1,  where  all  pij  >  0  and  s,  =  y pXJ  <  1. 

i= i 

The  last  condition  is  that  the  row  sums  are  less  than  or  equal  to  1.  If  all  row  sums  equal  1,  the 
matrix  is  said  to  be  stochastic;  if  all  row  sums  are  less  than  1,  the  matrix  is  said  to  be  properly 
substochastic. 

The  row  sums  have  a  significant  effec-  on  the  eigenvalues  of  P,  as  can  be  seen  from  the  following. 
Let  A  be  the  dominant  eigenvalue  of  P,  the  eigenvalue  of  maximum  modulus.  It  is  known  that 
A  <  1;  A  =  1  if  P  is  stochastic.  Suppose  x  =  (ii,...,x„)  is  the  left  eigenvector  belonging  to  A, 
normalized  so  that  x*  =  1.  Note  that  such  a  normalization  is  possible  since,  for  the  dominant 
eigenvalue,  all  X{  >  0.  Then  xP  =  Ax,  and  taking  row  sums  of  both  sides,  £  JNS.  =  A.  Therefore, 

I 

n 

/*  =  1-A  =  yx,(1-s,). 

1=1 


12 


Only  the  terms  for  which  s,  <  1  appear  in  the  sum.  If  P  is  indecomposable,  all  z,  >  0  and  we  see 
again  the  fact  that  if  P  is  indecomposable  and  properly  substochastic,  then  the  dominant  eigej  value 
A  <  1. 

Such  matrices  can  arise  from  finite  difference  approximations  of  Laplace’s  equation.  (As  such, 
applicability  to  a  wide  range  of  physical  problems  is  seen.)  Consider  a  square  grid  in  two  dimensions 
where  the  nodal  points  are  denoted  by  (i,j),  where  i,j  are  positive  integers.  If  u  is  a  function  on 
this  grid  Au  at  the  point  (i,j)  is  approximated  by 

^2  [ui+  l,j  +  ui-l,j  +  ti.J+l  +  ~  4 Uij), 

where  h  denotes  the  grid  size.  Setting  h  =  1,  one  sees  that  -Ac  =  I  -  A,  where  I  is  the  identity 
matrix  and  A  is  a  stochastic  matrix  whose  non-zero  entries  are  usually  1/4;  the  order  n  of  both 
matrices  is  the  number  of  grid  points.  Suppose  A  is  modified  so  that  each  boundary  grid  point  is 
an  absorbing  point.  Then  A  takes  the  following  form: 


h 
R  P 


Here  1^  is  the  identity  matrix  of  order  d,d  being  the  number  of  boundary  grid  points,  and  P  is 
a  substochastic  matrix  of  order  n  —  d.  In  fact,  P  must  be  properly  substochastic;  otherwise,  the 
boundary  would  not  be  connected  to  our  region.  The  number  of  rows  dx  of  P  which  have  sum  <  1 
can  now  be  interpreted  as  the  number  of  grid  points  directly  adjacent  to  the  boundary.  We  should 
expect  d\  <  d  with  near  equality.  Ideally,  equality  should  be  an  indication  of  smoothness  of  the 
boundary,  for  example,  the  absence  of  corners. 

It  will  be  noticed  that  other  algebraic  properties  of  P  correspond  to  geometric  properties  of  the 
region  which  is  modeled  by  the  grid.  For  example,  failure  of  the  region  to  be  connected  corresponds 
to  a  matrix  P  which  is  completely  decomposable,  that  is,  by  renumbering  indices  we  obtain  a  matrix 
of  the  form, 


Pa  0 
0  P22 


Indecomposable  matrices,  those  matrices  which  can  be  put  in  the  form 


Pit  P12 
0  P22. 


where  Pj2  ^  0,  simply  do  not  occur  in  this  context.  Additional  properties  of  reversible  matrices, 
and  matrix  reductions  have  been  derived  [16]. 

We  can  apply,  for  example,  the  foregoing  io  the  solution  of  the  following  dynamical  problem: 
let  us  consider  one  of  the  mplest  examples  in  mechanics,  the  harmonic  oscillator.  Let  x(t)  be  a 
function  for  t  >  0  which  is  a  solution  of  the  differential  equation 

x  +  u2x  =  0, 

where  x(0)  =  xq,  i( 0)  =  0,  so  that  x(t)  =  xocosu>t,  x(t)  =  -wzosin  uit.  If  w  is  held  constant, 
(x,-i)  describes  an  ellipse  in  phase  space. 

Since  the  orbit  is  closed,  we  can  approximate  the  system  by  a  finite  transition  probability 
matrix  as  follows.  Divide  the  curve  into  n  divisions  and  let  a,;  denote  the  probability  the  particle 
lands  in  the  jth  interval  in  a  unit  time  supposing  that  it  was  in  the  interval  just  before  the 
transition.  If  the  intervals  are  numbered  consecutively  1,2,3, ...,n,  so  that  after  n  we  regain  1, 
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then  the  matrix  A  =  (a^)"J=1  should  have  the  property  that  the  matrix  be  the  same  regardless 
of  which  interval  is  numbered  1.  Thus,  we  are  led  to  conclude  that  ai+  1<J-  =  a^j- 1.  Here  the 
convention  is  adopted  that  the  subscripts  are  taken  mod  n;  e.g.,  021  =  ai„.  Thus,  if  the  first  row 
is  given  by  au  =  Pi- 1 ,  for  i  =  1, . . . ,  n,  all  other  rows  are  given  by  atJ  =  pj_,  (again  with  the  mod 
n  convention  for  subscripts).  If  T  is  the  nth  order  matrix  defined  by 

0  10  0 
0  0  1  0  ••• 

0  0  0  1 


00  •••  1 
10  0  0 

then  A  is  given  by  A  =  52™=$  Pi T*.  The  eigenvalues  of  T  are  seen  to  be  Xj  =  exp(2xtj/n),  where 
j  =  0,  l,...,n  —  1  and  i  =  y/—l-  Indeed,  if  x  =  (xo,xi,. . .  ,x„-i)T  is  a  right  eigenvector  for  T, 
then  Xk  =  exp(2nijk/n).  These  must  still  be  the  eigenvectors  for  A,  but  now  the  eigenvalues  of  A 
are 

n  — 1 

A;  =  £  Pk  exp  (2,r‘Jfc/ n)< 

k=Q 

which  displays  each  Xj  as  a  mean  of  nth  roots  of  unity.  From  this  it  is  clear  that  two  cases  arise: 

1.  All  pk  except  one  are  equal  to  0.  In  this  case  A  is  periodic,  i.e.,  Am  =  I  for  some  m  which  is 
a  divisor  of  n.  Thus,  if  v0  is  an  initial  probability  distribution,  it  will  recur  after  m  time  steps: 
no  =  u0Am. 

2.  At  least  two  pk  >  0.  In  this  case  /imm_o0  Am  exists,  and  in  fact  we  must  have  u0Am  = 

n_1(l,  1,...,1)  for  any  arbitrary  initial  probability  distribution  v0.  In  other  words,  as  time 
goes  to  infinity  we  must  approach  a  state  where  there  is  equal  probability  of  the  particle  being 
in  any  of  the  n  subintervals. 

Note  that  case  1  is  in  fact  the  deterministic  case.  In  the  above  formulation  the  only  relevant 
numerical  data  is  the  probability  distribution  (po,...  ,Pn-i)  which  can  be  directly  calculated  from 
the  probability  distribution  for  the  natural  frequency  (or  spring  constant)  of  the  system.  This  can 
be  said  to  govern  the  speed  at  which  the  particle  goes  around  the  orbit:  T*  is  the  matrix  which 
corresponds  to  a  particle  moving  through  k  subintervals  in  a  unit  time.  Nothing  is  said  regarding 
the  radius  of  the  orbit  (actually  one  should  more  strictly  speak  of  the  semi-major  axis  since  the 
orbit  is  an  ellipse).  The  radius  of  the  orbit  is  a  function  of  the  energy,  which  is  constant  if  no 
external  force  is  applied.  Suppose  now  there  is  such  a  force  applied.  Then  we  may  suppose  its 
effect  is  to  alter  the  radius  of  the  orbit.  Since  the  loading  probability  may  be  assumed  independent 
of  the  structural  properties  of  the  system,  we  may  solve  this  problem  independently  of  the  former 
problem;  i.e.,  the  radial  distribution  is  independent  of  the  tangential  distribution. 

An  example  of  a  stochastic  matrix  of  order  n  +  2  for  such  a  transition  of  orbit  in  unit  time  is 
given  by 
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The  meaning  of  Q  is  as  follows:  there  exist  n  adjacent  orbits  and  two  absorption  states  such 
that  a  particle  has  equal  probability  of  transferring  to  one  of  its  two  adjacent  orbits  or  to  an 
absorption  level,  from  which  there  is  no  return.  Then  A  =  1  is  an  eigenvalue  of  Q  with  multiplicity 
2;  corresponding  left  eigenvectors  can  be  taken  as  (1,0, ...,0)  and  (0,...,0,1).  The  remaining 
eigenvalues  of  Q  are  the  eigenvalues  for  the  portion  of  Q  which  is  marked  off  in  the  dashed  box. 
Applying  to  P,  the  same  method  employed  earlier,  we  see  that  the  principal  right  eigenvector  for 
P  must  be  of  the  form,  a  =  (aj,...  ,an)T,  where  a*  =  sin(A:x/(n  -f  1)),  where  the  corresponding 
eigenvalue  is  given  by  A  =  cos  (x/(n  +  1)).  In  fact,  the  totality  of  eigenvalues  and  right  eigenvectors 
for  P  is  given  by  replacing  x/(n  +  1)  by  jx/(n  4-  1),  for  j  =  1,2, . . .  ,n.  Thus  if  Ai  is  the  dominant 
eigenvalue,  — Aj  is  also  an  eigenvalue,  and  actually,  all  eigenvalues  occur  with  their  negative  values. 
This  could  also  have  been  seen  from  the  fact  that  the  trace  of  P  equals  0.  Because  of  this,  it 
is  not  possible  to  use  our  previous  theory  to  associate  a  probabilistic  meaning  to  the  eigenvector 
associated  to  Aj. 

But  this  shortcoming  can  be  remedied  very  easily  by  considering  P2  instead  of  P.  P2  has 
positive  eigenvalues  but,  except  for  the  possible  A  =  0,  all  occur  with  multiplicity  2.  However,  P  is 
completely  decomposable  and  it  suffices  to  work  with  only  one  of  the  “parts”  of  P.  The  physical 
interpretation  illuminates  the  foregoing  algebraic  assertion:  in  two  time  steps,  the  system  must  lie 
in  an  even  numbered  state,  if  it  started  out  in  an  even  numbered  state,  with  the  same  if  we  replace 
“even”  by  “odd”.  Thus,  if  we  double  our  time  step  and  consider  only  the  even  numbered  states  we 
obtain  a  transition  matrix,  which  we  again  call  P: 


0  1  1 

U  4  2J 

This  new  P  should  have  an  order  approximately  half  that  of  the  previous  P,  but  again  we  may 
suppose  the  order  to  be  n.  The  dominant  eigenvalue  is  now  given  by 
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The  corresponding  left  (or  right)  eigenvector  is  given  by  the  components  a*  =  sin(fcx/[n  +  l]),k  = 

1,2 ,...,«. 

If  the  absorption  states  are  removed  and  n  permitted  to  go  to  infinity,  it  will  then  be  found 
that  the  limiting  distribution  approaches  that  of  a  Gaussian  or  normal  distribution. 


MARKOV  CHAIN  TRANSITION  MATRICES  as 
CONSTITUTIVE  MODELS  of  SOIL  DYNAMIC  BEHAVIOR 

To  apply  a  Markov  chain  model,  we  recall  that  the  transition  probabilities  {Xn}  with  state  space 
{0,1,2,...},  best  exhibited  in  the  form  of  a  matrix 
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must  be  derivable.  The  elements  of  a  transition  probability  matrix  P(m,n )  satisfy  the  conditions 


Pj,k(m,n)  >  0  for  all  j,  k 
2>.*(m,»)  =  1  for  all  j. 

k 


One  may  derive  the  recursive  relation: 


p(n)  =  p(0)P(0,n), 


where 


p(n)  =  lpo(n)  pi(n)...pj(n)...], 
Pj(n)  =  P[Xn=j}. 


It  follows  that  the  probability  law  of  a  Markov  chain  {An}  is  completely  determined  once  one 
knows  the  transition  probability  matrices,  and  the  unconditional  probability  vector  p(0)  at  time  0. 
In  the  case  of  a  homogeneous  (time-invariant)  Markov  chain  {Xn},  let 

P(n)  =  {pj,k(n)j,  P  =  {Pi,*} 

denote  respectively  the  n-step  and  the  one-step  transition  probability  matrices.  From  the  above 
equations,  it  is  observed  that 


P{n)  =  Pn, 
p(n)  =  p(0)Pn. 

Consequently,  the  probability  law  of  a  homogeneous  Markov  chain  is  completely  determined  once 
the  one-step  transition  probability  matrix  P  =  {Pj,k}  is  known,  and  the  unconditional  probability 
vector  p(0)  =  (pj(0)}  at  time  0. 

Markov  theory  is  a  conceptual  framework  for  modeling  the  propagation  in  time  and  space  of 
system  and  environmental  uncertainties.  As  a  theoretical  construct  it  is  extensive  [6,14,23].  When 
it  is  used  to  gain  physical  understanding  about  a  specific  dynamical  system,  it  is  necessary  that 
certain  parameters  particular  to  that  system  be  incorporated  into  its  mathematical  structure.  These 
parameters  will  represent  the  stochastic  material,  geometric,  and  environmental  characteristics  of 
the  system  being  studied.  The  transition  probability  matrix  is  the  vehicle  by  which  such  information 
is  incorporated  into  the  Markov  model. 

Others  have  taken  the  Markov  chain  approach  to  physical  modeling  [5,19,  for  example].  How¬ 
ever,  the  text  below,  tc  the  best  knowledge  of  this  author,  is  not  available  elsewhere.  An  Appendix 
summarizes  some  of  the  mportant  definitions  from  Markov  theory. 

Definition  of  States 

Implicit  in  the  above  discussion  is  the  fact  that  the  Markov  model  operates  on  a  system  which 
has  been  partitioned  into  states.  Inherent  in  the  definition  of  the  state  space  for  a  problem  is  an 
understanding  of,  even  if  only  in  general  terms,  what  are  the  likely  system  response  patterns  to 
the  environments  to  which  it  will  be  exposed.  The  state  space  is  mutually  exclusive  and  complete; 
all  possible  states  can  represented  in  a  unique  fashion.  It  is  desirable  to  devise  a  procedure  for 
state-space  specification  that  is  unbiased,  meaning  that  it  does  not  impose  the  modeler’s  a  priori 
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expectations  on  the  types  of  behavior  the  model  can  predict.  Ideally,  one  would  aim  for  a  model 
mechanism  which  would  generate  and  augment  the  state  space  using  initial  and  current  state 
information. 

The  process  of  state  space  definition  will  begin  with  the  determination  of  behavior  regimes 
that  are  of  particular  importance  or  interest.  Just  as  with  a  finite  element  model,  where  more 
elements  are  specified  in  regions  of  complexity  or  possible  critical  behavior,  the  definition  of  a  state 
space  exhibits  the  analyst’s  expectations.  Thus,  if  one  is  defining  the  state  space  for  a  material 
that  is  to  undergo  significant  strain  deformation  beyond  the  yield  strain,  ey,  one  might  define  a 
relatively  small  number  of  states  in  the  elastic  region  and  a  larger  number  in  the  yielded  region. 
Uncertainty  regarding  the  numerical  value  of  cy  will  require  the  definition  of  several  yield  states 
with  corresponding  transition  probabilities  according  to  a  probability  distribution  law  derived  from 
data.  Similarly,  uncertainty  about  the  elastic  loading  modulus  may  motivate  the  definition  of 
several  loading  moduli. 

Theoretically,  the  state  definition  procedure  is  not  a  trivial  task  since  the  analyst  must  antici¬ 
pate  behavior  states  that  may  not  have  been  previously  encountered,  and  do  this  without  hopelessly 
exceeding  analytical  and/or  computational  abilities.  It  is  also  not  an  isolated  task,  but  one  which 
must  be  carried  out  with  knowledge  about  the  ease  of  estimating  the  transition  probabilities.  It 
is  only  useful  to  define  states  which  can  be  observed  experimentally,  and  for  which  transition 
probabilities  can  be  estimated. 

In  particular  to  the  procedure  proposed  herein,  a  relatively  straightforward  definition  of  a  state 
space  and  derivation  of  transition  probabilities  based  on  experimental  data  is  demonstrated.  The 
following  example  displays  how  test  data  can  be  used  to  help  specify  the  state  space  for  a  process. 
Figure  1  depicts  the  results  of  a  series  of  10  dynamic  uniaxial  strain  tests  for  undisturbed  specimens 
of  CARES-Dry  sand  [7].  Differences  in  response  can  primarily  be  attributed  to  (small)  differences 
in  dry  density.  It  is  noted  that  the  stress-strain  curves  are  geometrically  similar,  but  spatially 
translated  along  the  strain  axis. 

A  possible  state  space  could  be  obtained  by  the  superposition  of  a  square/rectangular  grid  on 
the  above  figure  of  experimental  data.  This  approach  is  described  below. 

Estimation  of  Transition  Probabilities  From  Data 

The  main  thrust  of  this  section  is  now  developed;  that  is,  given  experimental  data  about  a  certain 
class  of  soil  behavior,  what  approach  can  be  used  to  derive  a  transition  probability  matrix  that  is 
a  reasonable  measure  of  the  evolutionary  properties  of  this  soil  media  class? 

Consider  the  generic,  parametric  xy  curves  of  Figure  2.  It  is  useful  to  be  able  to  derive  the 
transition  probability  matrix  for  the  process  which  results  in  this  behavior,  primarily  in  the  theory 
that  such  behavior  is  in  some  useful  sense  representative  of  near-media  behavior.  The  frequency 
interpretation  of  probability  will  be  used  here.  This  approach  assumes  that,  given  multiple  re¬ 
alizations  of  a  process  (as  one  would  have  due  to  a  sequence  of  experiments),  as  the  number  of 
realizations  becomes  large,  the  ratio  of  the  number  of  a  specific  realization  to  the  total  number  of 
realizations  is  a  reasonable  measure  of  the  probability  of  the  specific  occurrence.  This  approach 
assumes  also  that  each  realization  is  equally  probable,  although  one  could  appropriately  weigh  more 
likely  outcomes. 

Thus,  if  one  is  interested  in  estimating  the,  here  time-invariant,  transition  probability  between 
state  t  to  state  j,  7r*j ,  one  could  use  the  relation 
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where  the  denominator  represents  all  transitions  from  i,  including  »  — *  j,  and  those  that  are 
absorbed  within  i.  Using  this  procedure,  for  the  generic  experimental  data  of  Figure  2,  the  following 
transition  probability  matrix  II  can  be  derived  using  a  graphical  procedure: 


0 

11 

12 

13 

14 

22 

23 

24 

32 

33 

34 

0 

/ 

3/3 

11 

2/3 

1/3 

12 

1/2 

1/2 

13 

2/3 

1/3 

14 

1/2 

1/2 

22 

1/2 

1/2 

11=  23 

1/3 

1/3 

1/3 

24 

2/3 

1/3 

32 

1/1 

33 

1/3 

1/3 

34 

2/3 

43 

1/1 

44 

\ 

1/1 

43  44 

\ 


1/3 

1/3 


/ 


where  the  left  column  label  represents  the  initiating  state  (“from”),  and  the  top  row  label  represents 
the  receiving  state  (“to”).  It  is  assumed  that  from  the  initial  condition  0,  all  transitions  are  to  11. 
Thus,  the  |  entry  under  column  11. 

As  an  example  of  the  calculations,  consider  the  transitions  from  state  (2,2):  there  are  2  paths 
in  (2,2),  one  goes  to  (3,2),  the  other  to  (2,3);  thus,  each  is  assigned  a  probability  of  |. 

An  implicit  assumption  that  is  made  above  is  that  the  time  increments  between  transitions 
are  of  such  a  duration  that  transitions  will  only  occur  between  adjacent  (bordering)  states.  This 
requirement  can  always  be  satisfied  by  appropriately  refining  the  state  space.  In  general,  each  state 
may  lead  to  a  transition  to  any  of  8  states  for  the  uniform  state  space  grid  of  Figure  2. 

The  matrices  considered  here  are  sparse.  For  most  experimental  data,  this  will  be  the  case. 
Note  the  banded  nature  of  the  non-zero  elements.  This  is  also  a  result  of  the  manner  in  which  the 
state  space  is  numbered.  It  is  emphasized  that 

1.  the  derivation  of  the  above  transition  matrix  does  not  take  into  account  whether  the  individual 
curves  are  in  a  “loading”  or  “unloading”  phase;  this  distinction  will  be  important  for  some 
applications,  and  is  discussed  in  the  next  section,  and 

2.  while  the  transition  probability  matrix  is  homogeneous ,  an  implicit  time-scale  exists  due  to 
the  assumption  that  during  one  time  increment  only  the  adjacent  states  may  be  entered. 


Model  Refinement  for  More  Complicated  Behavior 


Here,  it  is  of  interest  to  model  behavior  such  as  the  reversal  of  trends  witnessed  in  stress-strain 
behavior  of  materials  subjected  to  loading  and  unloading  cycles.  The  procedure  of  the  previous 
section  is  modified  by  augmenting  the  state  space  such  that  “loading”  behavior  is  distinguished 
from  “unloading”  behavior.  Thus,  for  example,  in  Figure  2  state  (3,3)  is  in  actuality  2  states: 
(3,3/)  and  (3,3u),  where  the  suffix  /  denotes  loading,  and  u  denotes  unloading. 

With  this  approach,  one  can  partition  the  transition  probability  matrix  into  4  sub-matrices 
for  the  problem  at  hand: 

/  u 
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where  A  represents  an  evolution  from  loading  to  continued  loading,  B  represents  an  evolution  from 
loading  to  unloading,  C  represents  an  evolution  from  unloading  to  loading,  and  D  represents  an 
evolution  from  unloading  to  continued  unloading. 

For  the  generic  curves  of  Figure  2,  matrix  II  is  18  rows  by  18  columns,  and  the  sub-matrices 
are  as  follows: 
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and  C  =  0  in  this  instance,  where  it  is  noted  that  these  sub-matrices  are  sub-stochastic,  that  is, 
their  row  sums  are  <1.  Of  course,  the  row  sums  of  transition  matrix  II  are  equal  to  1;  II  is  a 
stochastic  matrix.  In  the  above,  all  rows/columns  exclusively  composed  of  0's,  except  for  the  “0” 
column,  have  been  removed. 

It  is  emphasized  that  the  above  transition  probability  matrix  is  only  representative  of  the 
“data  curves”  of  Figure  2,  and  is  a  product  of  the  state  space  assumed  in  that  figure.  This  means 
that  one  will  be  able  to  reproduce  these  curves,  using  the  derived  transition  probability  matrix,  by 
applying  the  theory  as  outlined  in  the  previous  section.  Thus,  one  would  obtain  the  following  state 
probabilities  at  each  increment  in  the  model: 
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p(  0)  =[1  00  0  0  0  0  0  000  000000  0] 

p(l)=  [0  1  0  0  0  0  0  0  000  000000  0] 

p{  2)=  [0  0  0.67  0.33  0  0  0  0  000  000000  0] 

p(3)=  [0  0  0  0.33  0.33  0.17  0  0.17  0  0  0  0  00000  0] 

p(4)=  [0  0  0  0  0  0.50  0.08  0.17  0.25  0  0  0  00000  0] 

p(ll)=  [0  00  0  0  0  0  0  00  0.72  0.28  0  0  0  0  0  0] 

p(16)  =  [0  00  0  0  0  0  0  00  0.99  0.01  0  0  0  0  0  0] 

Such  state  information  would  be  useful  in  a  comprehensive  computational  tool  where  the 
evolutionary  behavior  of  soil  samples  is  required.  It  is  conceivable  that  tables  of  such  information 
can  be  developed  for  different  samples  under  different  environmental  conditions.  As  such,  a  look-up 
procedure  can  be  established  in  the  main  program. 

The  hope,  and  this  is  the  subject  of  further  work,  is  that  the  given  data  is  also  representative  of 
process  behavior  in  some  “near  field”,  the  definition  of  which  will  be  dependent  on  the  correlation 
structure  of  the  media.  To  assess  this  property,  it  will  be  necessary  to  obtain  similar  data  for  field 
points  near  the  original  point,  and  to  establish  procedures  for  estimating  such  correlations.  Given 
the  results  of  the  present  work,  and  of  the  necessary  extensions,  it  is  believed  that  a  framework 
for  modeling  the  evolutionary  behavior  of  a  media  such  as  soil  is  possible  using  this  phenomeno¬ 
logical  point  of  view.  In  a  general  formulation,  each  state  may  be  decomposed  into  “loading”  and 
“unloading”  states,  or  whatever  behavior  class  is  most  appropriate. 

This  approach  to  estimating  the  transition  probability  matrix  can  be  used  in  problems  where 
not  much  data  is  available,  but  one  can  draw  on  engineering  judgment  to  set  bounds  on  the  behavior 
and  loading.  In  this  instance,  one  may  draw  upper  and  lower  bound  a  -  e  curves  and  “spread”  the 
probabilities  evenly  (or  with  some  weighting)  throughout  the  state  space.  Such  a  transition  matrix 
may  provide  a  rough  measure  of  the  evolutionary  soil  behavior. 

The  state  space  definition,  a  rectangular  grid  above,  is  purely  a  construct  of  the  analyst,  and 
cylindrical,  spherical,  or  other  geometries  may  be  more  useful  for  other  types  of  problems.  The 
transition  matrix  may,  of  course,  be  used  to  represent  a  3-dimensional  problem.  The  key  is  always 
an  adequate  data  set,  or  other  perhaps  subjective  information,  for  estimating  probabilities. 

Applying  this  graphical  procedure  to  the  experimental  CARES-Dry  sand  data  is  straightfor¬ 
ward,  but,  due  to  the  10  samples,  one  may  superimpose  a  refined  grid  on  the  sample  paths,  which 
are  overlapping  and  “close”.  Transition  probabilities  for  this  application  will  be  in  multiples  of 
The  state  space  grid  need  not  be  uniform,  and  may  be  more  refined  in  regions  of  critical  behavior. 
Again,  the  transition  matrix  will  be  banded  and  sparse.  If  these  sample  paths  are  characteristic  of 
the  underlying  soil  media,  then  the  derived  transition  probability  matrix  will  be  representative  of 
the  soil  behavior. 

For  behavior  which  is  more  complicated  than  that  presented  above,  such  as  cyclical,  hysteretic 
behavior  of  soils,  one  may  extend  the  state  space  definition  such  that  each  complete  loading¬ 
unloading  cycle  is  represented  by  a  state-space  in  two  dimensions,  and  each  subsequent  loading¬ 
unloading  cycle  extends  the  state  space  in  the  third  dimension.  Thus,  a  three  dimensional  grid  will 
be  representative  of  soil  undergoing  many  cycle  hysteresis.  In  such  a  manner,  so-called  “memory- 
dependent”  behavior  can  be  accounted  for  in  the  model.  The  key  to  the  success  of  such  an  approach 
is  the  availability  of  sufficient  data.  But  then,  this  is  the  key  to  any  model  purporting  to  be 
representative  of  some  aspect  of  the  physical  world. 
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In  soil  dynamics  application,  the  analyst  would  use  the  Markov  chain  state  transition  model  to 
estimate  the  “location”  of  the  soil  in  its  state  space.  This  information,  i.e.,  p( 0),. . .  ,p(16)  above, 
would  be  fed  back  to  the  dynamics  computer  code  whenever  such  information  is  needed  to  update 
the  system  dynamics.  Some  thought  would  be  needed  as  to  whether  the  most  probable  state  is 
used,  or  some  average  of  possible  states  be  used,  to  define  the  evolution  of  the  soil  in  its  state 
space.  Use  of  a  particular  transition  matrix  presupposes  that  a  similar  soil,  and  loading,  exist 
for  the  problem  at  hand.  It  is  intended  that  the  transition  matrix  derivation  be  automated;  given 
experiment -generated  data  points,  an  algorithm  can  be  used  to  generate  the  transition  probabilities 
for  specific  problems  in  real  time. 


CONCLUDING  REMARKS  &  SUGGESTIONS  for  ADDITIONAL  WORK 

We  have  learned  much  during  the  course  of  our  work  here.  Certain  aspects  of  our  work,  had  we 
been  able  to  retrace  our  steps,  would  have  received  less  emphasis,  and  other  asperte  more  emphasis. 
But,  of  course,  that  is  a  decision  that  can  only  be  made  a  fortiori.  Where  we  have  made  progress:  a 
generic  framework  for  modeling  complex  media  such  as  soils  is  identified  and  explored;  an  approach 
to  extract  probabilistic  information  from  experimental  data  is  proposed  and  demonstrated  using 
what  we  believe  is  the  most  comprehensive  soil  data  available.  Where  it  would  have  been  interesting 
to  have  reached  a  more  advanced  stage:  more  progress  would  have  been  desirable  in  translating 
the  mathematical  understanding  of  Markov  transition  matrices  into  a  physical  understanding  of 
the  constitutive/dynamic  behavior  of  the  soil. 

A  continuation  of  this  work  would  require  the  following  components.  The  first  is  the  exten¬ 
sion  of  some  of  the  theoretical  lines  drawn  in  categorizing  stochastic  matrices,  their  convergence 
properties,  and  the  linking  of  matrix  classes  with  dynamic  behavior.  In  addition,  the  derivation  of 
transition  matrices  which  represent  more  complicated  behavior,  such  as  hysteresis,  would  demon¬ 
strate  the  validity  of  the  graphical  methodology  presented  in  general.  As  part  of  this,  testing  the 
robustness  of  this  methodology  would  be  an  enlightening  activity.  Automating  the  procedure  so 
that  transition  matrices  can  be  generated  as  data  is  compiled  in  an  experiment  could  prove  very 
useful  in  linking  mathematical  model  to  data. 

One  issue  of  prime  concern,  and  one  which  becomes  paramount  to  those  who  must  work  with 
complicated  media  such  as  soil,  is  the  question  of  the  modeling/predictive  limits  of  analytical 
methods.  That  is,  as  good  as  our  mathematical  models  and  experimental  techniques  become,  there 
is  an  inherent  uncertainty  in  the  medium  that  is  of  such  significance  that  the  results  of  any  analysis, 
while  rigorous,  may  not  provide  information  which  will  be  useful  in  a  practical  (predictive)  sense 
(except  perhaps  over  extremely  short  time  durations). 

The  problem  alluded  to  above  is  a  function  of  the  “scales”  of  the  physical  problem  under 
study.  For  example,  a  site  of  1,000  cubic  feet  may  be  reasonably  modeled,  given  general  information 
regarding  the  site  environment,  and  some  minimal  tests  within  that  10  foot  cube.  This  is  because 
any  discontinuities  will  have  a  reasonable  chance  of  being  identified,  and  that  other  “continuous 
changes”  will  not  likely  be  significant  due  to  the  size  of  the  length  scale,  that  is,  the  length  scale  is 
small  relative  to  any  discontinuities. 

If  one  now  considers  a  cube  with  sides  of  length  1,000  foot  (lxlO9  cubic  feet),  a  different 
class  of  problem  arises.  This  is  because  discontinuities  can  be  completely  contained  within  the 
volume,  and  may  be  missed  by  tests.  The  significance  of  the  site  environment  is  less  important. 
In  addition,  continuous  changes  in  properties  now  will  work  over  a  larger  length  scale.  If  one  adds 
the  fact  that  all  variations  in  the  medium  have  random  components  that  are  easily  on  the  order  of 
the  mean  value,  then  it  can  be  concluded  that  any  known  model  of  the  medium  will  not  likely  be 
representative  of  even  adjacent  sites. 
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This  problem  of  “scales”  of  physical  behavior  is  encountered  in  other  disciplines,  such  as 
turbulence  in  oceanic  and  atmospheric  processes,  where  parameter  variations  over  the  scales  of 
the  prc'  !^m  are  significant,  and  therefore,  predictions  are  limited.  la  fact,  small  scale  turbulence 
results  m  large  scale  response.  Similar  difficulties  exist  for  the  problem  at  hand. 

How  does  one  transcend  this  problem?  These  concluding  comments  will  not  be  a  source  for 
such  resolution,  rather,  a  few  brief  observations.  Is  the  answer  with  microstructural  modeling? 
We  are  not  qualified  to  address  this  technical  approach.  However,  if  a  model  is  microstructural 
instead  of  phenomenological,  then  the  above  problem  of  scales  will  tend  to  become  more  intricate. 
For,  even  if  at  a  microscale  an  accurate  mathematical  model  is  derived,  this  model  must  be  able 
to  predict  gross  behavior  at  the  larger  scales  of  interest  in  application.  To  do  this,  the  generalized 
forces  and  pressures  that  a  microstructural  model  predicts  at  microscale  must  “accumulate”  in  some 
sense  to  the  physical  forces  and  pressures  measurable  in  the  laboratory  and  field.  To  do  this,  the 
microstructural  model  must  be  valid  at  many  orders  of  magnitude  of  behavior,  from  the  microscopic 
to  the  aggregate.  Otherwise,  the  microstructural  model  will  be  only  a  formal  expression. 

Considering  the  difficulties  of  models  which  attempt  predictions  at  much  smaller  behavior 
ranges  and  scales,  it  will  be  interesting  to  discover  how  microstructural  models  are  approaching  the 
difficulties  briefly  discussed  above. 
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Appendix.  FRACTAL  MODELING  OF  SOIL  MEDIA 


The  basic  purpose  of  this  section  is  to  investigate  whether  soil  microstructure  characteristics  can  be 
modeled  by  fractal  sets.  The  purely  geometrical  part  of  this  work  involves  the  setting  up  of  a  model 
that  simulates  the  soil  structure  in  appearance.  It  is  necessary,  however,  to  go  beyond  this  and  to 
devise  a  model  that  will  not  merely  resemble  the  soil  in  appearance  but  also  in  mechanical  behavior. 
By  this  is  meant  the  deformation  under  loading,  including  the  failure  mechanism  phenomenon. 

In  the  following,  a  procedure  is  summarized  for  calculating  the  deformation  under  loading  of 
a  fractal  set  model.  While  the  treatment  is  necessarily  brief,  there  is  sufficient  detail  devoted  to 
particular  simple  examples  to  demonstrate  the  feasibility  of  a  program  that  has  the  potential  of 
initiating  a  subject  that  can  become  as  voluminous  as  continuum  mechanics. 

It  is  the  usual  practice  in  studying  the  deformation  of  structural  materials  to  model  these 
materials  as  continua.  By  a  continuum  we  mean  a  material  whose  mass  and  elastic  (or  inelastic) 
parameters  vary  continuously  through  the  medium.  In  short,  no  voids  of  discontinuities  are  permit¬ 
ted  in  such  a  model,  except  a  manageably  small  finite  number  of  these.  This  permits  the  powerful 
methods  of  calculus  to  be  employed,  resulting  in  a  system  of  differential  equations,  the  solution  of 
which  lead"  to  a  complete  description  of  the  deformation  of  the  material  under  an  arbitrary  loading. 

Microscopic  examination  shows  that  no  material  can  be  described  as  a  continuum.  Thus,  a 
continuum  model,  if  it  works  at  all,  works  when  the  deformations  vary  continuously  when  averaged 
over  sufficiently  large  volumes.  This  can  lead  to  some  inadequacy  in  dealing  with  important  inelastic 
phenomena,  such  as  failure.  It  might  be  said  then  that  continuum  models  are  in  their  predominance 
today  due  to  mathematical  convenience  rather  than  to  physical  observation. 

Lately  the  fractal  set  has  gained  attention  in  that  many  structures  that  are  found  in  nature 
are  more  readily  modeled  as  fractal  sets  rather  than  continua.  In  this  work,  we  intend  to  show  how 
such  a  fractal  set  model  can  lead  to  a  theory  of  deformation  which  parallels  that  of  a  continuum 
material  model. 

We  will  forgo  generality  and  confine  ourselves  to  the  simplest  possible  example  that  displays 
the  basic  ideas.  This  will  enable  us  to  avoid  the  formidable  technicalities  which  surround  the  theory 
of  fractals.  Thus,  only  an  intuitive  notion  of  a  fractal  set  will  suffice  to  understand  our  approach. 

We  now  give  a  complete  solution  to  one  of  the  simplest  problems  imaginable:  the  fractal  analog 
of  the  one-dimensional  stress-strain  test.  Imagine  a  layered  medium,  which  is  composed  of  rigid 
layers  interspersed  with  elastic  layers.  To  convert  the  problem  into  one  of  fractal  mechanics,  we 
assume  the  medium  is  of  unit  depth  and  that  the  rigid  layers  occupy  precisely  the  locations  of  the 
middle-third  intervals. 

To  make  this  fractal  model  more  plausible,  one  may  consider  it  as  the  limit  of  a  succession  of 
mass-spring  models.  This  is  shown  in  Figure  3,  where  the  first  two  stages  are  shown.  In  the  first 
stage,  shown  in  Figure  3a,  the  total  mass  m  is  concentrated  in  a  single  block  of  length  1/3,  while 
two  springs  k  of  equal  length  make  up  the  rest  of  the  model.  Upon  further  resolution  shown  in 
Figure  3b,  the  springs  further  resolve  into  4  springs  of  length  1/32 ,  and  so  on.  In  the  limit,  the 
springs  will  occupy  no  length  at  all!  In  contrast,  consider  the  equivalent  continuum  model,  which  is 
conceived  as  being  the  limit  of  mass-spring  models,  as  shown  in  Figure  4.  In  the  limit  the  springs 
occupy  the  total  length  of  the  model,  while  the  masses  occupy  zero  length.  Viewed  in  this  manner, 
the  fractal  model  is  no  more  bizarre  than  the  continuum  model. 
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The  equations  for  the  continuum  approximation  illustrated  in  Figure  4b  are  given  as  follows: 

m 

+  «(2uj  -  u2)  =  Fj 
—  u2  +  k(—ui  +  2u2  +  U3)  =  F2 

771 

—  ii3  4-  fc(— u2  +  2U3)  =  F3  +  kA  (1) 

Let  us  consider  the  static  case  with  F\  =  F2  =  F3  =  0.  Then  the  system  becomes: 

2uj  -  u2  -  0 
~ux  +  2u2  +  u3  =  0 

-ti2  +  2u3  =  A  (2) 

for  which  the  solution  is  given  by  uj  =  ^,u2  =  ^,«3  =  |A.  If  we  rewrite  this  solution  as 
u/c  =  kA/2 2  for  k  =  1,2,3,  then  it  will  be  apparent  that  at  the  nth  stage,  when  there  are  2"  -  1 
masses,  the  solution  will  be  given  by  ti*  =  fcA/2"  for  k  =  1,2,..  .,(2n  -  1).  Since  the  masses  are 
equidistant,  tne  solution  approaches  the  linear  function  u(x)  =  kA. 

When  the  same  approach  is  applied  to  the  fractal  approximations  given  in  Figure  3b,  we 
obtain  exactly  the  same  equations  given  in  (2);  (note,  however  that  the  masses  in  the  dynamic 
version  of  (1)  would  be  slightly  different).  Thus,  the  solution  would  be  exactly  the  same.  But  the 
interpretation  of  u*  is  quite  different;  in  the  limit  the  solution  is  given  by  u(x)  =  /( x)A,  where  / 
is  known  as  the  Devil’s  Staircase  function  [17). 

As  a  matter  of  fact,  had  we  known  nothing  of  the  Devil’s  Staircase,  we  could  have  used  this 
example  to  define  it.  Pushing  this  line  of  thought  further,  we  can  define  other  fractal  functions. 
Refer  again  to  the  continuum  model  in  Figure  4.  Now  let  A  =  0  and  all  F)  =  F.  This  case  models 
a  rod  acted  upon  by  its  own  weight.  Proceeding  as  before,  but  omitting  details,  when  n  — ►  00,  we 
find  the  expected  solution,  that  u(x)  is  proportional  to  the  quadratic  function  x(l  -x).  If  we  follow 
the  same  procedure  in  the  case  of  the  fractal  model  shown  in  Figure  3,  we  will  obtain  a  fractal 
function  which  is  the  analog  of  the  quadratic  function.  In  like  manner  other  fractal  functions  can 
be  defined. 

The  first  thing  to  note  about  the  examples  above  is  that  both  the  continuum  and  fractal  models 
can  be  approximated  by  systems  having  a  finite  number  of  degrees  of  freedom,  and  when  they  are 
so  approximated,  the  resulting  equations  are  identical.  It  is  in  the  interpretation  of  the  solution  of 
these  equations  that  the  difference  lies. 

Consider  the  expression,  u*_i  -  2u*  +  Ui+\.  For  the  continuum  model,  this  expression  becomes 
proportional  to  the  second  derivative  of  u  at  some  point  on  the  model,  as  the  degree  of  approxi¬ 
mation  goes  to  infinity.  This  is  so  because  u;  represents  the  motion  at  a  point,  and  the  relevant 
points  are  equidistant. 

But  in  the  fractal  model,  represents  the  deflection  of  an  entire  interval,  and  no  such  inter¬ 
pretation  can  be  ascribed  to  u<_  1  -  2uj  +  u,+j.  Indeed,  the  limit  function  u  is  not  differentiable, 
and  it  is  meaningless  to  speak  of  a  fractal  model  being  governed  by  a  differential  equation. 

The  second  thing  to  note  is  that  for  the  fractal  model  solution,  the  deformation  function  u(z) 
is  given  by  a  fractal  function,  i.e.,  a  function  that  varies  only  on  a  fractal  set.  These  functions  may 
appear  strange  at  first,  in  that  they  are  nondifferentiable;  nevertheless,  they  are  quite  amenable 
when  they  lie  under  an  integral  sign. 

This  means  that  the  finite  element  method  may  be  used  in  conjunction  with  a  fractal  model 
much  the  same  way  as  with  continuum  models.  One  need  only  recall  how  the  finite  element  method 
works.  Over  some  simple  geometric  element,  say  a  triangle  or  tetrahedron,  a  restricted  class  of 
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deformations  is  allowed.  This  class  of  deformations  will  consist  of  polynomials,  the  coefficients 
of  which  are  linear  functions  of  the  deformations  of  the  vertices  of  the  elements.  For  a  fractal 
model,  one  simply  replaces  the  polynomials  by  suitable  fractal  functions.  In  other  words,  the  finite 
element  method  can  be  employed  for  fractal  models  in  much  the  same  way  as  for  continuum  models; 
polynomial  functions  are  replaced  by  fractal  functions. 

The  method  presented  above  is  very  general.  It  remains  to  elaborate  these  in  two  and  three 
dimensions.  This  will  involve  the  determination  of  fractal  functions  and  their  integration  over 
various  finite  elements.  Comparisons  with  continuum  models  will  be  useful.  These  ideas  can,  in 
principle,  be  generalized  to  more  complex  problems. 


APPENDIX.  THE  MARKOV  METHOD  [23] 


In  classical  physics,  a  basic  role  is  played  by  the  fundamental  principle  of  scientific  determinism: 
from  the  state  of  a  physical  system  at  the  time  to,  one  may  deduce  its  state  at  a  later  instant  t.  As 
a  consequence  of  this  principle  one  obtains  a  basic  method  of  analyzing  physical  systems:  the  state 
of  a  physical  system  at  a  given  time  <2  may  be  deduced  from  a  knowledge  of  its  state  at  an  earlier 
time  and  does  not  depend  on  the  history  of  the  system  before  time  tj.  An  example,  differential 
equations. 

For  physical  systems  which  obey  probabilistic  laws  rather  than  deterministic  laws,  one  may 
enunciate  an  analogous  principle:  the  probability  that  the  physical  system  will  be  in  a  given  state 
at  a  given  time  may  be  deduced  from  a  knowledge  of  its  state  at  any  earlier  time  t\ .  and  does 
not  depend  on  the  history  of  the  system  before  time  t\.  Stochastic  processes  which  represent 
observations  on  physical  systems  satisfying  this  condition  are  called  Markov  processes. 

A  special  kind  of  Markov  process  is  a  Markov  chain;  it  may  be  defined  as  a  stochastic  process 
with  an  evolution  which  may  be  treated  as  a  series  of  transitions  between  the  states  of  the  system. 
These  states  have  the  property  that  the  probability  law  of  the  future  evolution  of  the  process,  once 
it  is  in  a  given  state,  depends  only  on  this  state  and  not  on  how  the  process  arrived  at  this  state. 
The  number  of  possible  states  is  either  finite  or  countably  infinite. 

A  discrete  parameter  (time  in  this  case)  stochastic  process  {X(t),t  =  0, 1,2, ...}  or  a  continuous 
parameter  stochastic  process  {X(t),t  >  0}  is  said  to  be  a  Markov  process  if,  for  any  set  of  n  time 
points,  ti  <  *2i  •  •  •  ,tn,  the  conditional  distribution  of  X{tn),  for  given  values  of  X(ti), . . .  ,X(tn_i), 
depends  only  on  X(fn_ j),  the  most  recent  known  value;  in  mathematical  notation,  for  any  real 
numbers  ii, . . .  ,i„ ,  representing  possible  states  of  the  system, 

F(X(tn)  <  Xn\X(tx)  =  *1, .  . .  ,  X(tn-l)  =  In  — l]  =  -P[A(t„)  <  XnlA^n-j)  =  X„_i]. 

This  equation  may  be  read  to  mean  that,  given  the  “present”  of  the  process,  the  “future”  is 
independent  of  its  “past.” 

Markov  processes  are  classified  according  to  (i)  the  nature  of  the  index  set  of  the  process 
(discrete  or  continuous  time),  and  (ii)  the  nature  of  the  state  space  of  the  process  (discrete  or 
continuous  state). 

A  real  number  x  is  said  to  be  a  possible  value,  or  state,  of  a  stochastic  process  X(t)  if  there 
exists  a  time  t  such  that  the  probability  P[x  -  h  <  X(t)  <  x  +  h]  is  positive  for  every  h  >  0.  The  set 
of  possible  values  of  a  stochastic  process  is  called  its  state  space.  The  state  space  is  called  discrete 
if  it  contains  a  finite  or  countably  infinite  number  of  states.  A  state  space  which  is  not  discrete  is 
called  continuous.  A  Markov  process  with  discrete  state  space  is  called  a  Markov  chain.  Integers 
{0,1,...}  are  used  to  represent  the  state  space  of  such  a  chain. 

A  Markov  process  is  described  by  a  transition  probability  function,  often  denoted  by 
P(x,to’,E,t)  or  P(E,t\x,to),  which  represents  the  conditional  probability  that  the  state  of  the 
system  will  at  time  t  belong  to  the  set  E,  given  that  at  time  to  <  t  the  system  is  in  state  x.  The 
Markov  process  is  said  to  have  stationary  transition  probabilities,  or  to  be  homogeneous  in 
time,  if  P(x,  to;  E,  t)  depends  on  f  and  to  only  through  the  difference  (t  -  to).  We  initially  consider 
homogeneous  in  time  models. 

The  theory  of  Markov  chains  is  initially  presented  with  the  assumption  that  the  transition 
probabilities  and  matrices  are  known.  This,  of  course,  is  never  the  case;  the  estimation  of  these 
transition  probabilities  is  usually  the  most  difficult  and  elusive  part  of  the  analysis. 
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In  order  to  specify  the  probability  law  of  a  discrete  parameter  Markov  chain  {A„}  it  suffices 
to  state  for  all  times  n  >  m  >  0,  and  states  j  and  k,  the  probability  mass  function 


Pj(n)  =  P[Xn  =  j] 


and  the  conditional  probability  mass  function 

Pj,k(m,n)  =  P[Xn  =  k\Xm  =  j]. 

The  function  Pj,fc(m,n)  is  called  the  transition  probability  function  of  the  Markov  chain.  The 
probability  law  of  a  Markov  chain  is  determined  by  pj(n)  and  pJ)*(m,n)  since  for  all  integers  q, 
and  any  q  time  points  ni  <  ni  <  ...  <  nq,  and  states  ki,...,kq, 

^t^ni  —  ki, ,  Xn 4  =  kq]  —  Pk !  (ni)pfc,,fc3(ni,n2)  ...  Pk(q-l),kq{nq~l  i  nq)- 

A  Markov  chain  is  said  to  be  homogeneous  (or  to  be  homogeneous  in  time  or  to  have  sta¬ 
tionary  transition  probabilities)  if  pJifc(m,n)  depends  only  on  the  difference  n  -  m.  We  then  call 

Pj.k(n)  =  P[Xn+t  =  k\Xt  =  j]  for  any  integer  t  >  0 

the  n-step  transition  probability  function  of  the  homogeneous  Markov  chain  {Xn}.  In  words, 
Pjtk(n)  is  the  conditional  probability  that  a  homogeneous  Markov  chain  now  in  state  j  will  move 
after  n  (time)  steps  to  state  k.  The  one-step  transition  probabilities  p>,jt(l)  are  usually  written 
simply  pjtk,  i.e., 

Pj.fc  =  Ppft+r  =  k\Xt  =  j)  for  any  integer  t  >  0. 

Similarly,  if  {A'(t),t  >  0}  is  a  continuous  parameter  Markov  chain,  then  to  specify  the  probability 
law  of  {X(t),f  >  0},  it  suffices  to  state  for  all  times  t  >  s  >  0,  and  states  j  and  k ,  the  probability 
mass  function 

Pk(t)  =  P[X(t)  =  k] 

and  the  conditional  probability  mass  function 

Pj,k(s,t)  =  P[X(0  =  *IX(s)=j). 

The  function  pjtk(s,t)  is  called  the  transition  probability  function  of  the  Markov  chain.  The 
Markov  chain  is  said  to  be  homogeneous  (or  to  have  stationary  transition  probabilities)  if  Pj,*(a,f) 
depends  only  on  the  difference  t  -  s.  We  then  call 

Pi,k(t)  =  P[X(t  +  u)  =  fc)A'(u)  =  j]  for  any  u>0 

the  n-step  transition  probability  function  of  the  Markov  chain  X(t),t  >  0. 

A  fundamental  relation  satisfied  by  the  transition  probability  function  of  a  Markov  chain  {A„} 
is  the  so-called  Chapman-Kolmogorov  equation:  for  any  times  n  >  u  >  m  >  0  and  states  j 
and  k, 

PjAm’n)  =  ]Tp>,.(m-u)p<,*(u>n). 

summing  over  all  states  t  of  the  Markov  chain.  This  equation  represents  all  possible  intermediate 
states  »  between  j  and  k. 
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•••  Pi,/fc(m,n) 
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•••  Pj,k(min) 

. 

PJ,K(m,n). 

The  transition  probabilities  of  a  Markov  chain  {Xn}  with  state  space  {0,1,2,...}  are,  best 
exhibited  in  the  form  of  a  matrix: 


P(m,n)  = 


Note  that  the  elements  of  a  transition  probability  matrix  P(m,n)  satisfy  the  conditions 

pJifc(m,n)  >  0  for  all  j,k 
y;  Pj,k(m,  n)  =  1  for  all  j. 

k 

The  Chapman-Kolmogorov  equations  for  all  times  n  >  u  >  m  >  0  may  be  written: 

P(m,n)  -  P(m,u)P(u,n). 

From  the  Chapman-Kolmogorov  equation,  one  may  derive  the  recursive  relation: 

p(n)  =  p(0)P(0,n), 

where 


P(n )  =  [Po(n)  Pi(n) . .  >Pj(n) . . .], 

Pj(n)  =  P[Xn  =  j]. 

It  follows  that  the  probability  law  of  a  Markov  chain  {X„}  is  completely  determined  once  one 
knows  the  transition  probability  matrices,  and  the  unconditional  probability  vector  p(0)  at  time  0. 
In  the  case  of  a  homogeneous  Markov  chain  {X„},  let 

P(n)  =  {Pj.it(n)},  P  =  {pj.fc} 

denote  respectively  the  n-step  and  the  one-step  transition  probability  matrices.  From  the  above 
equations,  it  is  observed  that 

P(n)  =  P", 
p(n)  =  p(0)Pn. 

Consequently,  the  probability  law  of  a  homogeneous  Markov  chain  is  completely  determined  once 
one  knows  the  one-step  transition  probability  matrix  P  =  {pj,k},  and  the  unconditional  probability 
vector  p( 0)  =  {pj(0)}  at  time  0. 

A  Markov  chain  is  said  to  be  a  finite  Markov  chain  with  K  states  if  the  number  of  possible 
values  of  the  random  variables  {Xn}  is  finite  and  equal  to  K.  The  transition  probabilities  p^,*  are 
then  non-zero  for  only  a  finite  number  of  values  of  j  and  k ,  and  the  transition  probability  matrix 
P  is  then  a  KxK  matrix. 

A  state  k  is  said  to  be  accessible  from  a  state  j  if,  for  some  integer  N  >  1,  Pj,k(N )  >  0.  Two 
states  j  and  k  are  said  to  communicate  if  j  is  accessible  from  k,  and  k  is  accessible  from  j. 

A  state  k  is  said  to  be  recurrent  if  the  probability  is  1  that  the  Markov  chain  will  eventually 
return  to  k ,  having  started  at  k.  A  state  k  is  said  to  be  non-recurrent  if  the  above  probability 
is  less  than  1. 

A  state  k  is  called  an  absorbing  state  if  pk,k  —  so  that  once  the  chain  visits  k  it  remains 
there  forever.  An  absorbing  state  is  clearly  recurrent. 
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Figure  1.  WES  Data  [7] 
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Figure  3.  Fractal  Model 
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Figure  4.  Continuum  Model 
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